arXiv:1505.02379vl [cond-mat.mtrl-sci] 10 May 2015 


On the accuracy of commonly used density functional approximations in determining 
the elastic constants of insulators and semiconductors 

M. RasandeiQ and M. A. Moram 

Department of Materials, Imperial College London, SW7 2AZ, London, United Kingdom 

(Dated: May 12, 2015) 

We have performed density functional calculations using a range of local and semi-local as well as 
hybrid density functional approximations of the structure and elastic constants of 18 semiconductors 
and insulators. We hnd that most of the approximations have a very small error in the lattice 
constants, of the order of 1%, while the error in the elastic constants and bulk modulus are much 
larger, at about 10%. In addition, we find that the error in the elastic constants, dj, are larger 
compared to the error in the bulk modulus. Depending on the functional and which error estimate 
that is being used, the difference in the error between the elastic constants and the bulk modulus 
can be rather large, about a factor of two. According to our study, the overall best performing 
density functional approximation for determining the structure and elastic properties is the PBEsol, 
closely followed by the two hybrid functionals PBEO and HSE, and the AM05 functional. 


I. INTRODUCTION 


Fundamental material parameters are needed to de¬ 
sign, characterize and simulate new devices effectively. 
For example, accurate elastic constants are essential for 
determining the composition of epitaxial films using X- 
ray diffraction^ for assessing the critical thicknesses for 
strain relaxation in device heterostructured 2 *and for mod¬ 
eling the behavior of dislocations in these materials.^ 
Elastic constants can also be used to obtain the geometry 
of the individual layers in- and out-of-plane in multilayer 
structures! 4 However, obtaining elastic constants exper¬ 
imentally is sometimes difficult. For the Ill-nitrides, for 
example, it is very difficult to determine the elastic con¬ 
stants accurately through experiment due to the inher¬ 
ently large measurement uncertainties and limitations in 
the sample quality. Therefore, theoretical elastic con¬ 
stants values are used routinely for these systems.®' Fur¬ 
thermore, for other s emi conductor systems such as many 
II-IV nitride systemtP^ theoretical elastic constants are 
absolutely necessary since experimental values are lack¬ 
ing. It is therefore necessary to have theoretical tools 
that are accurate and efficient with the additional re¬ 
quirement to have a predictive power. Density func¬ 
tional theory (DFT) is such a tool and it is now estab¬ 
lished that the crystal structure, lattice dynamics and 
electronic structure of almost any element or compound 
can be accurately treated by this theory. Neglecting the 
technical and numerical aspects of the implementation, 
the accuracy in a DF calculation is determined by the 
exchange-correlation (XC) energy functional. The first 
approximation for the XC energy functional that was de¬ 
veloped is the local density approximation (LDA) which 
assumes that the XC energy of the electron density in a 
material is identical to the XC energy of the free elec¬ 
tron gas with the same electron density. Even though 
the LDA is a very simple approximation, it has been 
found to give remarkably good results for many mate¬ 
rial properties. However, this level of accuracy is not 
sufficient for many of the applications that have been 


mentioned above. Many approximations have however, 
been proposed in order to improve on the LDA, starting 
with the jjeneralisedgradient approximations (eg. PBEp 
PW91 Ml revPBEpJ BFB eP P BEsoP 1 ) or other semi¬ 
local functionals (eg. A MO fP3 0®3), to the more elaborate 
meta-GGAs (eg. TPSSP® 1 3) and hybrid density func¬ 
tionals (eg. PBElPand HSEP®). 

The aim with the present study is to establish the level 
of accuracy of many of the previously mentioned XC en¬ 
ergy functional approximations, namely the LDA, PBE, 
AM05, PBEsol, RPBE, TPSS, PBEO and HSE approxi¬ 
mations, in determining the structure and especially the 
elastic properties of semiconductors and insulators. Pre¬ 
vious comparative studies for solids have m ostly f ocused 
on the descriptio n of the latt ice constant d 15 l 17 l 21 l tH anc j 
the bulk modulu^^i 7 ^2212il f or both metallic and semi¬ 
conducting/insulating systems. We will discuss both the 
accuracy in the lattice constants and the bulk modulus. 
Furthermore, we will discuss the accuracy in determining 
the single crystal elastic constants, c^. The bulk mod¬ 
ulus is an important physical property of a material. It 
measures the response of a material under hydrostatic 
compression. However, it contains less information than 
the knowledge of the elastic constants. It is therefore 
important to assess the accuracy of DF calculations in 
determining the elastic constants in relation to the ac¬ 
curacy in determining the structural properties and the 
bulk modulus. Especially, for the applications listed pre¬ 
viously where the individual elastic constants are abso¬ 
lutely essential. 

The 18 systems that have been investigated in this 
work include elemental semiconductors (C, Si and Ge), 
zinc-blende semiconductors (BN, BP, GaP, GaAs, InP, 
InAs, InSb, and SiC) and insulators in the rock-salt struc¬ 
ture (LiF, LiCl, NaF and MgO). Additionally, we have 
included two systems in the fluorite structure (CaF 2 and 
Mg 2 Si) as well as the small band gap skutterudite CoSb 3 . 
These latter systems have been included to provide more 
structural complexity compared to the traditional zinc- 
blende and rock-salt systems. Especially, CoSb3 has a 
large unit cell containing 16 atoms, where large voids are 
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present in the structure, see for example Ref. [27! and ref¬ 
erence therein. 

The paper is outlined as follows: In Section [XT] we will 
give a brief theoretical background on the DF approxi¬ 
mations which have been used in this study. For more 
details, we refer to the original publications. In Section 
|III| we present the details of our calculations and in Sec¬ 
tion [TV] we will present our results. Finally, in Section |V| 
we will summarise our findings and present the conclu¬ 
sions. 

II. THEORETICAL BACKGROUND 


LDA gives reasonably reliable geometries for solids and 
for elastic constants the typical error has in previous 
studies been found to be in the order of 5-109PSI de¬ 
pending on the system. However, it fails badly for the 
atomisation energies of molecules and solids. 

The generalized gradient approximation (GGA) con¬ 
stitute a family of approximations where the XC energy 
can be expressed as 

Eg GA [n] = J ri(r)e GGA (n(r), Vn(r))d 3 r, 

= J n(r)e^ A (r s {r))F xc (r s {r),s{r))d 3 r , (4) 


Within the Kohn-Sham approach to density functional 
theory, the total energy of a system of electrons is given 
by (all equations are given in atomic units) 


Etot = T S + v ext (r)n(r)d 3 r 


" (r) ” (r VdV + E„, (1) 


r — r' 


where T s is the kinetic energy of a system of noninter¬ 
acting electrons and the following terms represent the 
electron-nucleus energy, electron-electron energy and E xc 
is the exchange-correlation (XC) energy which is un¬ 
known and has to be approximated as will be discussed 
below. This term may be divided further into an ex¬ 
change and a correlation term: E xc = E x + E c . 


A. Local and semi-local functionals 

For a local or semi-local density functionals, the XC 
energy functional can, in general, be described by 

E xc [n-f, ri|] 

= /dV„(r )e „(„ T . % ,V„ f ,V % , r„ n ), (2) 

where n a , Vn a and r a are the electron density, gradi¬ 
ent of the density and the kinetic energy density of spin 
a =f or respectively. However, we will from now on 
neglect the spin index. Within the local density approx¬ 
imation (LDA) the XC energy only depends on the elec¬ 
tron density and can be expressed as 

E X c A [n } = J n(r)e^ A (n(r))d 3 r, (3) 

where the XC energy per unit volume e x ® A is a func¬ 
tion of the electron density n(r) and is chosen to be 
identical to the XC energy of the uniform electron gas 
with the same density. The exchange part is given by 
4 DA = — (3/4)(3/7r) 1 / 3 n 4 / 3 . The correlation energy is 
obtained through a fit^ of accurate quantum Monte 
Carlo calculations for the uniform electron gasP The 


where F xc (r s ,s) = F x (s) + F c (r s ,s) is the enhancement 
factor, r 3 = 3/(47ra) is the Wigner-Seitz radius and 
s = |Vn|/[2(37r 2 ) 1 / 3 n 4 / 3 ] is the reduced density gradi¬ 
ent. In the literature two types of GGA functionals can 
be found: (i) the empirical functionals, whose parame¬ 
ters were determined by fitting to experimental or first 
principles data, and (ii) t he parameter free functionals, 
e.g. PBEP and PW91®^, whose parameters were deter¬ 
mined in order to satisfy mathematical relations which 
are known to hold for the exact functional. Within this 
context it should be mentioned that the parameter free 
functionals contain arbitrary choices such as the analyt¬ 
ical form chosen to represent F xc or the choice of con¬ 
straints to be satisfied. Here we will only discuss param¬ 
eter free functionals. 

The PBE functional® was designed to satisfy several 
conditions obeyed by the exact functional, e.g. the cor¬ 
rect uniform electron gas limit (i.e. LDA is recovered 
when s = 0), the Lieb-Oxford bouncP ( E x > E xc > 
— 1.679 f ?r 4 / 3 d 3 r) and the LDA linear response. The en¬ 
hancement factor for exchange is given by 


F™{s) 


1 + K — 


K 


( 5 ) 


where n = 0.804 and /r = 0.21951. One advantage of 
the PBE functional is that it performs well for finite 
and infinite systems and therefore rather diverse in its 
applications. For solids, the PBE has a tendency of 
underbinding, i.e. yielding too large binding distances 
and too small formation energies, and several functionals 
have been designed to improve on the performance of the 
PBE for solids. Here we will briefly mention two, namely 
PBEsoP and RPBEP. 

The PBEsoP was designed to be more accurate than 
PBE for solids and surfaces. Within a GGA there is a 
choice to be made between the accuracy in e.g. lattice pa¬ 
rameters and surface energies in comparison to accurate 
atomisation energies. The PBE, as has been mentioned, 
performs well for both solids and molecules, however, an 
improved functionality in one area of applications, e.g. 
for solids, will worsen the performance for atoms and 
molecules. Within the PBEsol the same analytic form as 
in the PBE is used, with some of the parameters changed 
to satisfy modified constraints in comparison to PBEP 
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TABLE I. Mean error (ME), mean absolute error (MAE), root mean square error (RMSE) and mean absolute relative error 
(MARE) for lattice constants, a , elastic constants, dj, averaged over all combinations of ij (i.e. ij = 11,12 and 44), and for 
the bulk modulus, B, for all functionals used in this work. The values shown in bold are the measures that show the closest 
agreement with the experimental values. 


a (A) eg (GPa )B (GPa) 


XC 

ME 

MAE 

RMSE 

MRE 

MARE 

ME 

MAE 

RMSE 

MRE 

MARE 

ME 

MAE 

RMSE 

MRE 

MARE 

LDA 

-0.047 

0.047 

0.061 

-0.9% 

0.9% 

7.1 

9.2 

15.9 

5.9% 

8.1% 

7.2 

7.9 

10.9 

8.0% 

8.8% 

PBE 

0.076 

0.076 

0.086 

1.4% 

1.4% 

-8.7 

11.6 

14.7 

-8.7% 

10.8% 

-8.9 

9.0 

11.1 

-8.9% 

8.9% 

AM05 

0.024 

0.033 

0.044 

0.5% 

0.6% 

-2.9 

7.8 

11.5 

-4.3% 

7.4% 

-3.1 

5.9 

7.0 

-4.1% 

6.3% 

PBEsol 

0.013 

0.025 

0.031 

0 . 3 % 

0 . 5 % 

- 1.5 

7.4 

11.6 

- 2 . 1 % 

6.8% 

- 0.9 

5.1 

6.5 

- 0 . 7 % 

5.3% 

RPBE 

0.135 

0.135 

0.149 

2.6% 

2.6% 

-15.6 

17.2 

21.8 

-15.2% 

16.1% 

-16.2 

16.2 

18.1 

-16.6% 

16.6% 

TPSS 

0.052 

0.052 

0.061 

1.0% 

1.0% 

-6.1 

10.1 

13.5 

-6.2% 

9.1% 

-6.6 

7.7 

9.2 

-5.9% 

7.8% 

PBEO 

0.023 

0.025 

0.033 

0.4% 

0 . 5 % 

6.6 

8.4 

16.3 

3.0% 

5.5% 

4.0 

4.5 

8.2 

2.4% 

3.2% 

HSE 

0.025 

0.030 

0.039 

0.5% 

0.6% 

5.9 

8.2 

16.0 

2 . 1 % 

5 . 3 % 

3.1 

4.1 

8.0 

1.2% 

2 . 7 % 


The PBEsol has indeed been found to yield improved lat¬ 
tice parameters and surface energies compared to PBE, 
however, the performance for atoms and molecules are 
worse than for the PBE. 

In order to achieve a XC functional with improved ad¬ 
sorption energies compared to LDA and PBE, Hammer 
et all 12 proposed a slightly adapted form of the exchange 
enhancement factor, F x , compared to the PBE, with 

F x rpbe (s) = 1 + k( 1 - e"^ 2/K ), (6) 

with an identical value for k and all correlation terms 
identical as in PBE. It was found that adsorption ener¬ 
gies of several atoms and molecules on metallic surfaces 
was greatly improved with the RPBE .^However, lattice 
constants are overestimated in RPBE and the bulk mod¬ 
ulus is underestimated significantly, which is related to 
the strong s dependence compared to the PBEP^ 

While the standard procedures for developing new 
GGA functionals are to fit parameters to experimental 
data and/or by satisfying universal mathematical condi¬ 
tions on the exact DF, the AM05 functional was the first 
functional to use a subsystem functional schemed The 
idea is to use separate functionals from different model 
systems for which the XC energy is known, namely the 
uniform electron gas and the surface jellium. A DF index, 
which depends on the reduced density gradient, is then 
used to locally determine the nature of the systemJ^l^ 
AM05 is therefore a systematic improvement over LDA, 
by the inclusion of terms that depend on the density gra¬ 
dient while maintaining the XC limit of the LDA. It has 
been found that the PBEsol and AM05 have very simi¬ 
lar performance for many systems, with the PBEsol per¬ 
forming slightly better than AMOSP^This was argued by 
Csonka et al 2A to be due to the very similar behaviour 
of the XC enhancement factor F xc (s , r s ) between PBEsol 
and AM05 in solids where s < 1 everywhere. For some 
solids with s max » 1 the difference between PBEsol 
and AM05 is greater.^ 

If the dependence of the electron density, gradient of 
the density and the kinetic energy density are all used 
in the construction of the XC energy functional as de¬ 
scribed in Eq. Q we have a so called meta-GCA XC 


energy functional. Compared to a GGA functional the 
inclusion of the kinetic energy adds extra flexibility in the 
construction of a XC energy functional. Such functionals 
has the potential to correctly treat effects that cannot 
be treated accurately in LDA and GGA. Several types 
of meta-GGA functionals have been presented. Here, we 
will use the TPSS meta-GGA functional,^ which satis¬ 
fies several exact constrains without using any empirical 
parameters and it h as bee n found to be reliable for both 
molecules and solids ^ 


B. Hybrid density functionals 

Hybrid DF are non-local theories where some of the 
exchange part of a standard XC functional, e.g. PBE, 
has been substituted by some amount of Hartree-Fock 
exchange energy. The two hybrid functi onals that have 
been used in this study, PBE(jL^ and HSEjMAO both have 
25% Hartee-Fock exchange energy substituting regular 
PBE exchange. In the PBEO approximation the XC en¬ 
ergy is expressed by 

£tc BE 0 - |^ HF + ^, PBE + £ c PBE , (7) 

where E FBF and E B BE are the exchange and correlation 
terms in the PBE approximation and E BF is the Hartree- 
Fock exchange en ergy. In the HSE approximation the XC 
energy is given bjPSjESl 

£ BSE = E X SE + E bse , (8) 

where 

E bse = aE BF,SB (u]) + (1 - a)E™ E - SR M 

+ js ™,l* ( w)| ( 9 ) 

where LR and SR denote the long range and short range 
parts of the exchange energy respectively, see Refs. [19 ! 
and [201 a is a mixing parameter governing the amount 
of the non-local Hartree-Fock exchange and ui is a screen¬ 
ing parameter that controls the spatial range over which 
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FIG. 1. (Color online) Calculated relative error in the lattice 
constants, a, RE = (a — oo)/ao, where ao is the experimental 
reference. The bars are the relative errors for (from left to 
right) the LDA, PBE, AM05, PBEsol, RPBE, TPSS, PBEO 
and HSE approximations. 


FIG. 3. (Color online) Calculated relative error in C 12 , 
RE = (ci 2 — Ci 2 )/ci 2 , where cj 2 is the experimental reference. 
The bars are the relative errors for (from left to right) the 
LDA, PBE, AM05, PBEsol, RPBE, TPSS, PBEO and HSE 
approximations. Note that there is no experimental C12 data 
for CoSb 3 . 



FIG. 2. (Color online) Calculated relative error in cn, 
RE = (cn — Cii)/c?r, where Cn is the experimental reference. 
The bars are the relative errors for (from left to right) the 
LDA, PBE, AM05, PBEsol, RPBE, TPSS, PBEO and HSE 
approximations. 



FIG. 4. (Color online) Calculated relative error in C44, 
RE = (C44 — C44)/C44, where C44 is the experimental reference. 
The bars are the relative errors for (from left to right) the 
LDA, PBE, AM05, PBEsol, RPBE, TPSS, PBEO and HSE 
approximations. 


III. DETAILS OF THE CALCULATIONS 


the non-local exchange part is important. The correla¬ 
tion energy in the HSE approximation, E ese , is taken 
to be identical to the PBE correlation energy as in the 
PBEO approximation. The amount of Hartree-Fock ex¬ 
change in the HSE is identical to the amount in PBEO, 
i.e. a = 1/4. The range separation parameter, w, on the 
other hand has to be determined by comparison with ex¬ 
perimental data. It has been found that u> = 0.2—0.3 A -1 
gives good results with regards to structural as well as 
electronic properties of mate rials, with w = 0.2 A -1 be¬ 
ing the optimum choice.^ 2 °l 32 l We have in this study 
used lo = 0.2 A” 1 . Note that for w = 0 Eq. ^ is equiv¬ 
alent to the PBEO and, in addition, Eq. (f8l) asymptot¬ 
ically reaches the PBE for w —> oo! 1H l 33 T Hybrid DF 
are more accurate than LDA and PBE for many prop¬ 
erties. Especially, it is possible to obtain reliable band 
gaps using hybrid theories. However, the non-local na¬ 
ture of these approximations make them computation¬ 
ally very expensive. In addition, it has been found that 
PBEsol and AM05 are just as accurate as hybrid approx¬ 
imations for structural properties and in determining the 
bulk modulus.^ 


Density functional calculations have been performed 
using the projector augmented wave (PAW) method^ as 
implem ented in the Vienna ab initio simulation package 
(VASP) ! 35136 1 The plane wave energy cut-off was set to 
800 eV. For the local and semi-local density functional 
approximations (LDA, PBE, AM05, RPBE and PBEsol) 
we have used P-centered k-point meshes with the small¬ 
est allowed spacing between k-points of 0.1 A -1 and for 
the much more computationally demanding hybrid den¬ 
sity functionals (PBEO and HSE) the spacing was set to 
0.4 A -1 . For the meta-GGA TPSS functional we have 
used 0.1 A -1 . The atomic positions and simulation cell 
shapes were relaxed until the Hellmann-Feynman forces 
acting on atoms were smaller than 0.001 eV /A. 

Apart from when using the LDA, the PAW core po¬ 
tentials used in the calculations were obtained using the 
PBE functional. For most of the approximations, this 
is the most obvious choice since the approximations are 
extensions or revisions of the PBE. The only exception 
is AM05, but it has been showiP 5 that the difference 
between using LDA and PBE PAW core potentials in 
an AM05 calculation is small. For most elements in 
our study the PAW core potentials were set up using 
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FIG. 5. (Color online) Calculated relative error in the bulk 
modulus, B, RE = (B — Bq)/Bo, where Bo is the experimental 
reference. The bars are the relative errors for (from left to 
right) the LDA, PBE, AM05, PBEsol, RPBE, TPSS, PBEO 
and HSE approximations. 


obvious core-valence partitions. The exceptions were 
Ca (3s 2 3p 6 4s 2 ), Ga (3d 10 4s 2 4p 1 ), Ge (3d 10 4s 2 4p 2 ), In 
(4d 10 5s 2 5p 1 ) and Na (2s 2 2p 6 3s 1 ), where semi-core states 
have been treated as valence states. 

The elastic constants were evaluated following Refs.f37 
and [331 where the stress is evaluated from the application 
of a strain to the system and the elastic constants are 
evaluated from Hooke’s law, er = Ce, where er is the 
stress tensor, e is the strain tensor and C is the elas¬ 
tic constants tensor. This approach in combination with 
semi-local XC functionals has been used to calculate the 
elastic constants of both zinc-blendc 39 and wurtzit^HI 41 ! 
group Ill-nitride alloys, producing accurate el astic con - 
stants with a relatively low computational costEZMlM2J 
In order to assess our calculated data, these have been 
compared to available experimental data. Lattice and 
elastic constants have, as much as possible, been ex¬ 
tracted from low temperature experiments, ideally ex¬ 
tracted from T = 0 K. However, this is not always pos¬ 
sible in which case experimental values taken at room 
temperature have been used. For the lattice constants, 
we have whenever possible used the same values as Matts- 
son et where the same policy was used. In the case 
of CoSb 3 we were unable to find any measured value for 
the C 12 elastic constant and our calculated values for C 12 
presented here should be regarded as predictions. 

The accuracy of the calculated results have been eval¬ 
uated by the mean error (ME), the mean absolute error 
(MAE), the root mean square error (RMSE), the mean 
relative error (MRE) and the mean absolute relative error 
(MARE). The ME is given by 

n 

ME = — ^2 (Yi - Y i ), (10) 

i—1 

where Yi and Y) are the calculated and experimental val¬ 
ues for the quantity of interest, for example the lattice, 
a, or elastic, c^, constants, respectively, and n is the 
number of systems. Furthermore, the MAE is given by 


the RMSE is given by 


RMSE = 


\ 


-XF-F 2 , 

n L ^ 1 1 


the MRE is given by 


MRE = - V 

n 


i=l 


- Yi) 
Yi 


and the MARE is given by 


1 " \Yi-Yi 

MARE = V J-A-- 

n ^ Yi 


( 12 ) 


(13) 


(14) 


It will be shown that the choice of the error estimate will 
affect the conclusions slightly. Note that functionals that 
consistently over- or underestimate either the lattice con¬ 
stants or the elastic constants will have |ME| = |MAE| 
and | MRE | = |MARE|. 


IV. RESULTS AND DISCUSSION 


In TableUwe show the evaluated error measures for the 
lattice constants, elastic constants and the bulk modulus 
for the 18 solids in our study. Note that in Table[I]the er¬ 
ror estimates in the three elastic constants are averaged, 
i.e. A Cij = (niAcn+?r 2 Aci 2 + n 3 Ac 44 )/(ni+n 2 + n 3 ), 


where Ac,j is any of the measures ME, MAE, MRE or 
MARE, and n t is the number of evaluated systems for the 
elastic constant c^. Furthermore, note that the RMSE 
needs a slightly different treatment, which is apparent 
from the form of the RMSE in Eq. (12). Experimental 
references and the calculated results are shown in Ta¬ 
bles [TTJ [TTTJ [TVj [V] and [Vi] In addition, Figs. [T] through [5] 
show the relative error in the lattice constants, the elastic 
constants (one figure for each of the Cn, C12 and c 44 ) as 
well as the bulk modulus for each of the 18 systems in 
our study. 


A. Lattice constants 


Experimental lattice constants are usually measured 
at room temperature and, furthermore, they also include 
zero-point phonon effects (ZPPE). Therefore, these ex¬ 
perimental values are not directly comparable with the 
results of ground-state density functional calculations 
since the theory are obtained at T = 0 K. For lattice 
constants the ZPPE manifest as a zero-point anharmonic 
expansion (ZPAE) of the lattice.^ This effect can be es¬ 
timated bjfel 


MAE = 


1 n 

n 


(ii) 


AoQ Xpt _ AU 0 expt _ 3 (B . k B & D 

a e X pt gyjxpt 16 (^1 1 BgVQ Xpt ’ 


(15) 
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TABLE II. Lattice constants, a, for systems with cubic symmetry. The experimental lattice constants without the ZPAE 
subtracted are given within parenthesis. Data shown in bold show the least deviation from experimental values. All units in A. 



ExptP 

LDA 

PBE 

AM05 

PBEsol 

RPBE 

TPSS 

PBE0 

HSE 

c 

3.545 (3.567) 

3.533 

3.570 

3.549 

3.553 

3.587 

3.568 

3.543 

3.543 

Si 

5.415 (5.430) 

5.403 

5.470 

5.437 

5.436 

5.501 

5.456 

5.435 

5.438 

Ge 

5.646 (5.658) 

5.627 

5.761 

5.678 

5.674 

5.815 

5.705 

5.676 

5.682 

BN 

3.597 (3.615) 

3.581 

3.624 

3.604 

3.606 

3.587 

3.620 

3.595 

3.595 

BP 

4.522 (4.538) 

4.492 

4.549 

4.519 

4.521 

4.575 

4.546 

4.522 

4.523 

GaP 

5.439 (5.451) 

5.395 

5.506 

5.441 

5.437 

5.556 

5.487 

5.458 

5.463 

GaAs 

5.637 (5.648) 

5.610 

5.751 

5.669 

5.661 

5.808 

5.714 

5.674 

5.680 

InP 

5.855 (5.869) 

5.827 

5.957 

5.885 

5.876 

6.015 

5.948 

5.893 

5.899 

InAs 

6.048 (6.058) 

6.027 

6.185 

6.099 

6.086 

6.254 

6.159 

6.096 

6.103 

InSb 

6.473 (6.479) 

6.450 

6.631 

6.537 

6.518 

6.714 

6.600 

6.534 

6.542 

SiC 

4.341 (4.358) 

4.332 

4.380 

4.356 

4.359 

4.402 

4.369 

4.348 

4.349 

LiF 

3.969 (4.010) 

3.906 

4.074 

4.054 

4.022 

4.157 

4.049 

4.021 

4.023 

LiCl 

5.071 (5.106) 

4.963 

5.154 

5.129 

5.077 

5.266 

5.123 

5.118 

5.124 

MgO 

4.186 (4.207) 

4.153 

4.242 

4.211 

4.206 

4.285 

4.223 

4.193 

4.193 

NaF 

4.577 (4.609) 

4.507 

4.686 

4.664 

4.628 

4.810 

4.621 

4.645 

4.648 

CaF 2 

5.443 (5.466) 

5.316 

5.495 

5.437 

5.407 

5.591 

5.469 

5.447 

5.452 

Mg 2 Si 

6.329 (6.347) 

6.258 

6.359 

6.326 

6.324 

6.411 

6.355 

6.324 

6.329 

CoShP 

9.041 (9.055) 

8.910 

9.105 

8.966 

8.977 

9.174 

9.050 

9.033 

8.998 

ME 

_ 

-0.045 

0.077 

0.025 

0.014 

0.136 

0.053 

0.025 

0.026 

MAE 

- 

0.047 

0.077 

0.036 

0.027 

0.136 

0.053 

0.028 

0.033 

RMSE 

- 

0.061 

0.089 

0.047 

0.034 

0.152 

0.064 

0.037 

0.042 

MRE 

- 

-0.9% 

1.4% 

0.5% 

0.3% 

2.6% 

1.0% 

0.4% 

0.5% 

MARE 

- 

0.9% 

1.5% 

0.7% 

0.5% 

2.6% 

1.0% 

0.5% 

0.6% 


and the ZPAE c orre cted experimental value is OQ Xpt — 
ActQ Xpt . In Eq. (15), UQ Xpt is the experimental volume 
per atom, B\ is the the pressure derivative of the bulk 
modulus at equilibrium, B 0 is the bulk modulus at equi¬ 
librium, and Qd is the Debye temperature. For this cor¬ 
rection we have, as done previously in the literaturej 17 l 25 l 
used experimental values for Bq and @£>, while we have 
used theoretical values for B\. In Ref. [Tfl and [55] the¬ 
oretical Bi values were obtained using the TPSS func¬ 
tional, however, for reasons that will become apparent 
below, we have chosen to use the PBEsol DF instead. 
The pressure derivative of the bulk modulus, B i, were 
obtained by fitting 14 to 25 energy points evenly dis¬ 
tributed around the equilibrium volume to the equation 
of state of Vinet et aZ“ We find that the ZPAE correc¬ 
tion to the lattice constants is on average 0.019 A with a 
maximum for LiF (Ado = 0.041 A) and a minimum for 
InSb (Aao = 0.006 A). The errors shown in Tables [i] and 
|n]are calculated using the corrected values for the lattice 
constants. In Table HH we show both the corrected and 
uncorrected experimental values. In Table VII we show 
the experimental Debye temperatures and the calculated 
Bi values used in evaluating the corrections. 

Overall, for any DF approximation the error in the lat¬ 
tice constants is less or equal to 1.4%, except for the case 
of RPBE where the error is 2.6%. It is clear from Table [I] 
that the PBEsol is the best performing DF approxima¬ 
tion regarding the lattice constants, irrespective of the 
error measure. The PBEsol is closely followed by the 
PBE0, HSE and AM05 DF approximations. All of these 


approximations have a MARE smaller than 0.6%. The 
LDA and PBE approximations have a MARE of about 
1% and together with the RPBE and TPSS approxima¬ 
tions are the functionals that either overestimate or un¬ 
derestimate the lattice constants compared to experiment 
for all systems. The LDA is the only approximation that 
yields too short lattice constants for all systems. We find, 
as was also discussed by Csonka et aZ.pS that if the ZPAE 
correction had not been applied this would have given a 
bias to systems that have a tendency of yielding too long 
lattice constants, i.e. making, e.g. the PBE and RPBE 
approximations more accurate than they are, while mak¬ 
ing the LDA less accurate. We also find that the LDA 
is more accurate than the PBE, which was also found 
by Csonka et a/P when analysing their smaller set of 
10 non-metal solids. For their set of 14 metals the PBE 
performed better with a MARE of 1.2% while the LDA 
had a MARE of 2.7%P 

In Fig. [l] we show the relative error for each individ¬ 
ual system for the different DF approximations. It can be 
clearly seen that the LDA underestimates the lattice con¬ 
stants for all systems. Furthermore, the PBE, RPBE and 
TPSS approximations overestimates the lattice constants 
of all systems. The other approximations, only underes¬ 
timate the lattice constant for a few systems, and then by 
a very small amount. We note that regarding the lattice 
constants some systems are more difficult to achieve good 
agreement between experiment and theory, e.g. SiC and 
MgO. LDA underbinds the most for the rock-salt struc¬ 
tured crystals as well as for the two fluorite structured 
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compounds. 


B. Elastic constants 

As can be seen in Tables [TJ |III[ |IV| and [V] the errors in 
the elastic constants are larger, compared to the error in 
the lattice constants, with MARE of the order of 10% or 
better depending on the functional. In general, the DF 
approximations underestimate the elastic constants. It is 
only the LDA that consistently overestimate the values 
for the elastic constants. The PBEsol, PBEO and HSE 
approximations overestimate cn but underestimate both 
C12 and C44. We also find that the MAE and RMSE for 
cn are rather large, see Table m so that the MAE and 
RMSE errors in cn is larger than the corresponding er¬ 
rors in both C12 and C44 for all DF approximations. In 
addition, both the MAE and RMSE result in a larger er¬ 
ror in C44 compared to C12, except for the MAE for LDA 
which gives that the error in C 12 is about the same size 
as the error in C44. However, the absolute error derived 
from the MAE or RMSE is not the best way to obtain 
a proper view of the performance of the DF approxima¬ 
tions, since the size of the values for cn, C12 and C44 for 
each system is very different; typically Cn is much larger 
than both C12 and C44, and C44 is larger than Ci2- If we 
instead focus on the MARE, we find that the error in cn 
is smaller than the error in C12 and the error in C44 is 
smaller than the error in C 12 . In general, for the MARE 
it can also be concluded that the error in Cu is larger 
than the error in C44, with the exception being the PBE, 
AM05 and PBEsol approximations. It is because of the 
large variations in the size of the errors in the different 
elastic constants we also evaluated the averaged errors as 
shown in Table |T| The averaged error estimates for the 
elastic constants shown in Table U is intended to show 
the error that is expected for any single elastic constant 
when using any of the DF approximations. 

If we are to determine which functional that performs 
the best, we note that it depends on which measure that 
is being used. In general, however, for cn we find that 
PBEsol is the best performing functional, followed by 
the two hybrid approximation and AM05. For C 12 , the 
best performer is PBEO followed by HSE, PBEsol and 
AM05. Interestingly, for C44 LDA has the smallest ME, 
MAE, RMSE and MARE, while the HSE has the smallest 
MARE. For the averaged errors in Table |T[ PBEsol has 
the smallest ME, MAE and MRE, AM05 has the smallest 
RMSE, while HSE has the smallest MRE and MARE. 
However, the differences between PBEsol, AM05 and the 
two hybrid functionals regarding the elastic constants are 
overall not very significant. 


C. Bulk modulus 

It is common practise to evaluate the bulk modu¬ 
lus when comparing the performance of different den¬ 


sity functional approximations!^- 17 ^ The traditional 
method for determining the bulk modulus, H, is to calcu¬ 
late the total energy, E, as a function of volume, V, and 
fit the volume dependence of the energy to an equation 
of state. The bulk modulus at equilibrium can thereafter 
be obtained either as one of the fitting parameters or by 
evaluating 


B = 


1 d 2 E 

vW 2 ' 


(16) 


It is also possible to calculate the bulk modulus in terms 
of the elastic constants; for a cubic system it can be eval¬ 
uated bj0H 


B = 


cn + 2ci2 


(17) 


Here we have calculated the bulk modulus according to 
Eq. 

We note that it is possible to perform a similar correc¬ 
tion to the experimental bulk modulus as was done for 
the lattice constants.® However, when considering the 
bulk modulus, as well as the elastic constants discussed 
in the previous section, it is necessary to point out that 
(i) the experimentally determined bulk modulus can have 
a measurement uncertainty in the order of 10%, and (ii) 
the deviation between theory and experiment, as will be 
discussed, are much larger than for the lattice constants. 
Depending on the system and which DF approximations 
that is used the error varies from a couple of percent to 
about 20%. The ZPAE correction is therefore of less im¬ 
portance for the bulk modulus and the elastic constants. 
We have therefore compared our theoretical data with 
experimental data without corrections. 

The overall best performing functional is PBEsol, with 
the smallest ME, RMSE and MRE, while HSE has the 
smallest MAE and MARE. It is interesting to compare 
the error in the bulk modulus with the error in the elastic 
constants. In general, as shown in Table [TJ the error in 
the elastic constants are larger than the error in the bulk 
modulus. The MAE and RMSE are both consistently 
larger for the elastic constants. The ME gives more or 
less similar error for both the elastic constants and the 
bulk modulus while there is no particular trend for the 
ME and MRE. The MARE is larger for the elastic con¬ 
stants, except for the LDA and RPBE approximations, 
where the error in the elastic constants is smaller than 
the corresponding error in the bulk modulus. However, 
if we look at the best performing functional, the error 
in the elastic constants are almost twice as large as the 
error in the bulk modulus, see Tables [TJ |III[ [TVj |y| and 

EH 

The error in the bulk modulus can be written in terms 
of the errors in Cn and C 12 by using Eq. 

AB = -Acn + ~ Aci 2 ■ (18) 

From the above expression, it is clear that for a particular 
system the error in the elastic constants cn and C 12 will 






TABLE III. Elastic constant cn for cubic systems. Data shown in bold show the least deviation from experimental values 
All values are given in GPa. 



ExptP 

LDA 

PBE 

AM05 

PBEsol 

RPBE 

TPSS 

PBEO 

HSE 

c 

1079 

1106 

1056 

1079 

1075 

1032 

1052 

1145 

1143 

Si 

167 

161 

153 

156 

156 

149 

157 

171 

170 

Ge 

129 

123 

105 

116 

116 

99 

115 

133 

131 

BN 

820 

826 

781 

796 

795 

1032 

786 

845 

845 

BP 

315 

358 

339 

346 

346 

330 

340 

368 

367 

GaP 

140 

142 

124 

133 

134 

117 

126 

143 

142 

GaAs 

120 

116 

98 

107 

109 

91 

103 

118 

117 

InP 

100 

100 

87 

94 

94 

81 

87 

103 

101 

InAs 

83 

85 

71 

77 

79 

65 

72 

87 

86 

InSb 

67 

67 

54 

60 

61 

50 

56 

68 

67 

SiC 

390 

404 

384 

390 

390 

375 

388 

418 

417 

LiF 

136 

166 

117 

119 

134 

94 

121 

132 

131 

LiCl 

59 

78 

53 

54 

63 

40 

59 

57 

56 

MgO 

306 

340 

279 

297 

304 

252 

296 

317 

317 

NaF 

97 

136 

100 

98 

109 

79 

118 

105 

103 

CaF 2 

164 

191 

159 

164 

171 

144 

160 

171 

170 

Mg 2 Si 

121 

125 

117 

120 

120 

113.5 

118 

128 

127 

CoSbcpl 

164 

219 

183 

206 

211 

170 

195 

203 

204 

ME 

_ 

15.8 

-11.0 

-2.4 

0.6 

-22.9 

-6.0 

14.2 

13.2 

MAE 

- 

17.7 

16.1 

10.8 

10.7 

25.3 

14.5 

15.0 

14.5 

RMSE 

- 

24.3 

18.5 

15.5 

15.9 

29.3 

17.3 

24.1 

23.7 

MRE 

- 

9.4% 

-7.2% 

-2.6% 

0.4% 

-14.6% 

-3.4% 

4.9% 

4.1% 

MARE 

- 

10.8% 

9.7% 

6.6% 

6.5% 

15.6% 

8.7% 

5.8% 

5.5% 


have different contributions to the error in the bulk mod¬ 
ulus. An error in C 12 will have twice the contribution to 
the error in the bulk modulus compared to Cn- However, 
only one third of the error in cn, and two thirds of the er¬ 
ror in C 12 , is transferred to the error in the bulk modulus. 
A large error in cn and a small error in C 12 will thereby 
results in an error in the bulk modulus that is smaller 
than the error in the elastic constants. The same is also 
true for the opposite case of a small error in Cn and a 
large error in C 12 . On the other hand, if Acn ss Aci 2 this 
leads to A B « Acn, i.e. the error in the bulk modulus is 
the same as the error in c\\. However, even for the cases 
were the errors in Cn and C 12 are of similar size, the error 
in the bulk modulus is much smaller for the best perform¬ 
ing functionals, such as PBEsol, AM05, HSE and PBEO. 
The reason for this lies in an overall cancellation of errors 
due to the overestimation of cn and underestimation of 
C 12 for the PBEsol, HSE and PBEO. 


V. SUMMARY AND CONCLUSIONS 

We have performed DF calculations of the lattice con¬ 
stants, elastic constants and bulk modulus for a set of 
18 semiconductors and insulators. We find, in agreement 
with previous studies, that the overall best performing 
functional is PBEsol, followed by the two hybrid approx¬ 
imations PBEO and HSE, and AM05. These functionals 
are distinct improvements over the LDA and PBE ap¬ 
proximations. It should be kept in mind that the PBEsol 


and AM05 are from a calculation point of view much 
more efficient than PBEO and HSE due to the very large 
computational cost of the hybrid functionals. If a reli¬ 
able description of features of the electronic structure, 
e.g. band gaps, is required it is necessary to use the hy¬ 
brid functionals or to solve for the electronic structure 
by means of the GW approximation on top of a standard 
DF calculation. If structural and elastic properties are of 
interest, the PBEsol and AM05 functionals are the better 
choice for optimal efficiency, especially for large systems. 
Interestingly, the LDA performs better than the PBE for 
the lattice constants, bulk modulus and most of the elas¬ 
tic constants. It is only for cn were the error is smaller 
for the PBE approximation. 

The errors in the lattice constants are generally very 
small, less than 1.4%. It is only RPBE which gives a 
larger error. The errors in the elastic constants and bulk 
modulus, on the other hand, are much larger compared 
to the error in the lattice constants, of about 10% or 
smaller. Furthermore, we find that the error in the elastic 
constants are larger than the error in the bulk modulus. 
If the best performing functional is compared for both 
elastic constants and the bulk modulus, we find that the 
error in the elastic constants are about twice as large as 
the error in the bulk modulus. This is due to an over¬ 
all cancellation of errors between the cn and C 12 elastic 
constants. 

Finally, we note that the large deviations obtained us¬ 
ing the RPBE should not discourage the use of this func¬ 
tional. It was designed for improving adsorption energies 
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TABLE IV. Elastic constant C12 for cubic systems. Data shown in bold show the least deviation from experimental values. 
All values are given in GPa. Note that there is no experime ntal val ue fo r the C12 elastic constant of CoSb 3 . Th e value of C12 
can be evaluated to 45 GPa by using cn and B from Tables III and VI respectively, and by the use of Eq. (171. 



ExptP 

LDA 

PBE 

AM05 

PBEsol 

RPBE 

TPSS 

PBE0 

HSE 

c 

124 

149 

125 

141 

141 

115 

122 

138 

138 

Si 

65 

65 

57 

61 

62 

53 

59 

64 

63 

Ge 

48 

48 

37 

42 

44 

33 

42 

46 

44 

BN 

190 

191 

168 

181 

182 

115 

173 

186 

185 

BP 

100 

85 

73 

80 

81 

68 

73 

79 

79 

GaAs 

53 

53 

41 

46 

49 

36 

44 

50 

49 

GaP 

62 

64 

52 

57 

60 

46 

54 

61 

61 

InP 

56 

57 

46 

50 

53 

41 

47 

55 

54 

InAs 

45 

48 

37 

42 

44 

33 

40 

46 

45 

InSb 

37 

37 

28 

31 

34 

24 

29 

35 

34 

SiC 

142 

143 

128 

138 

138 

121 

132 

141 

141 

LiF 

47 

52 

49 

45 

47 

45 

46 

50 

50 

LiCl 

20 

23 

22 

20 

21 

21 

20 

23 

23 

MgO 

95 

98 

92 

91 

92 

88 

90 

99 

100 

NaF 

24 

24 

22 

21 

22 

21 

21 

23 

23 

CaF 2 

44 

61 

41 

43 

48 

32 

41 

46 

45 

Mg 2 Si 

22 

27 

23 

25 

25 

21 

24 

25 

25 

CoSb 3 

- 

47 

37 

42 

46 

33 

38 

33 

27 

ME 

_ 

3.0 

-7.8 

-3.5 

-1.8 

-12.9 

-6.9 

-0.5 

-1.0 

MAE 

- 

5.0 

8.5 

5.9 

4.8 

13.0 

7.1 

4.0 

4.3 

RMSE 

- 

8.6 

11.2 

7.9 

7.0 

15.8 

9.6 

6.5 

6.8 

MRE 

- 

6.3% 

-10.7% 

-5.9% 

-2.2% 

-18.9% 

-9.6% 

0.3% 

-0.9% 

MARE 

- 

8.6% 

13.0% 

9.0% 

7.2% 

19.4% 

10.5% 

6.3% 

6.7% 


compared to the PBE and for such applications it is very 
successful. 12 
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HSE 

c 

578 

594 

560 

582 

577 

547 

560 

613 

611 

Si 

80 

76 

74 

74 

74 

74 

77 

82 

82 
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TABLE VI. Bulk modulus for cubic systems. Data shown in bold show the least deviation from experimental values. All 
values are given in GPa. 



ExptP 

LDA 

PBE 

AM05 

PBEsol 

RPBE 

TPSS 

PBE0 

HSE 

c 

442 

468 

436 

454 

452 

421 

432 

473 

473 

Si 

99 

97 

89 

93 

94 

85 

92 

100 

99 

Ge 

75 

73 

60 

67 

68 

55 

66 

75 

73 

BN 

400 

403 

373 

386 

387 

420 

377 

405 

405 

BP 

172 

176 

162 

169 

169 

155 

162 

175 

175 

GaP 

88 

90 

76 

82 

85 

70 

78 

89 

88 

GaAs 

75 

74 

60 

66 

69 

54 

64 

73 

71 

InP 

71 

71 

59 

65 

67 

54 

61 

71 

70 

InAs 

58 

60 

49 

53 

56 

44 

50 

60 

59 

InSb 

47 

47 

37 

41 

43 

33 

38 

46 

45 

SiC 

225 

230 

213 

222 

222 

206 

217 

234 

233 

LiF 

77 

90 

71 

70 

76 

62 

71 

78 

77 

LiCl 

33 

41 

32 

31 

35 

27 

33 

34 

34 

MgO 

165 

179 

154 

160 

163 

143 

159 

172 

172 

NaF 

49 

61 

48 

47 

51 

40 

54 

51 

50 

CaF 2 

84 

104 

81 

84 

89 

69 

81 

88 

87 

Mg 2 Si 

55 

60 

55 

56 

57 

21 

55 

60 

59 

CoSb 3 

85 

104 

85 

97 

101 

79 

90 

90 

86 

ME 

_ 

7.2 

-8.9 

-3.1 

-0.9 

-16.2 

- 6.6 

4.0 

3.1 

MAE 

- 

7.9 

9.0 

5.9 

5.1 

16.2 

7.7 

4.5 

4.1 

RMSE 

- 

10.9 

11.1 

7.0 

6.5 

18.1 

9.2 

8.2 

8.0 

MRE 

- 

8 . 0 % 

-8.9% 

-4.1% 

-0.7% 

-16.6 

-5.9% 

2.4% 

1.2% 

MARE 

- 

8 . 8 % 

8.9% 

6.3% 

5.3% 

16.6% 

7.8% 

3.2% 

2.7% 


TABLE VII. Experimental Debye temperatures, Qn, and 
calculated first pressure derivative of the Bulk modulus, B i, 
for all systems in the present study. 



0d (Kjsi 


C 

2250 

3.68 

Si 

645 

4.27 

Ge 

373 

4.70 

BN 

1700 

3.67 

BP 

985 

3.78 

GaP 

445 

4.48 

GaAs 

344 

4.62 

InP 

425 

4.67 

InAs 

280 

4.78 

InSb 

160 

4.87 

SiC 

1200 

3.90 

LiF 

732 

4.33 

LiCl 

429 

4.40 

MgO 

945 

4.15 

NaF 

430 

4.72 

CaF 2 

510 

4.56 

Mg 2 Si 

417 

4.04 

CoSb 3 

307 

4.80 












